Spontaneous Patterning of Confined Granular Rods 
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Vertically vibrated rod-shaped granular materials confined to quasi-2D containers self organize 
into distinct patterns. We find, consistent with theory and simulation, a density dependent isotropic- 
nematic transition. Along the walls, rods interact sterically to form a wetting layer. For high rod 
densities, complex patterns emerge as a result of competition between bulk and boundary alignment. 
A continuum elastic energy accounting for nematic distortion and local wall anchoring reproduces 
the structures seen experimentally. 
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From nanometer sized molecules forming liquid crys- 
tals to timber floating down a river, rod shaped materials 
pervade our everyday existence. Onsager first demon- 
strated that hard core steric interactions suffice to order 
rods in thermal equilibrium 0. Above a critical con- 
centration, orientational entropy is sacrificed for gains in 
translational entropy; therefore, crowding alone induces a 
temperature-independent isotropic-nematic (I-N) transi- 
tion 2]. In this entropically driven transition, the kinetic 
motion of the rods serves as a mechanism for adequately 
sampling phase space. While thermal excitation mixes 
microscopic particles, macroscopic rods require an exter- 
nally applied energy for randomization. For both ver- 
tical and horizontal vibrations, a density dependent ne- 
matic to smectic-like transition was previously observed 
for granular rods Q. Such driven dissipative systems, 
however, are far from equilibrium and need not evolve 
to states of maximal configurational entropy. For exam- 
ple, prior studies utilizing macroscopic rods focused on 
jammed or other metastable states jHQ- 

Here, we report on finite-sized driven dissipative granu- 
lar systems of rods in steady-state that share many prop- 
erties of lyotropic liquid crystals at equilibrium. We find 
that short ranged hard core interactions determine the 
arrangement of the moving rods for all densities. A phase 
transition from a disordered isotropic state to a nematic- 
like state occurs as a function of rod density and rod as- 
pect (length-to-diameter) ratio, L/D. By using excluded 
volume scaling predicted by mean field theories for liq- 
uid crystals, our data give a single normalized transition 
density for all rod L/D and densities tested. 

However, excluded volume interactions in the bulk can- 
not explain all rod behavior because confining bound- 
aries also strongly influence rod ordering. For relatively 
low rod densities, experimental results for rod alignment 
and density profiles with respect to the wall strikingly 
resemble theoretical predictions |^], indicating that rods 
interact via hard core steric repulsions with the surround- 
ing boundaries. As the rods become dense, rod ordering 
in the bulk nematic competes with rod ordering along 
the surrounding walls. This frustration creates distinct 
yet easily reproduced patterns. We account for observed 



structures by modeling the system as a liquid crystal 
undergoing the elastic deformations of splay and bend 
subject to a simple wall interaction Moreover, we 

show that the patterns depend on the relative magnitude 
of wall and bulk energies. 

Rod-shaped granular materials confined to quasi-2D 
containers were vertically vibrated. Stainless steel wire, 
0.08 cm diameter, was cut into rods oi L/D ^ 20, 40, and 
60. Containers were either circular (radius i? = 15 cm or 
R — 7 cm) or square (29 cm diagonal) . The chamber 
height to rod diameter ratio was 20, which restricted the 
rods to lie mainly fiat. Due to the quasi-2D nature of 
these experiments, rod densities were scaled to the area 
of the container bottom, A, so that area fraction is = 
NLD/A, where N is the number of rods. A monolayer 
of rods lying flat and fully covering the container bottom 
corresponds to (p — 1. 

The rods were shaken with a sinusoidal acceleration of 
50 Hz and a peak acceleration of 4 times gravity. The 
magnetic field of the shaker was minimized by distancing 
the rod chamber from the shaker. Electrostatic forces 
were reduced by using aluminum chamber components 
and applying an anti-static coating to non-metallic sur- 
faces. The apparatus was carefully leveled. To randomize 
initial conditions, rods were dropped into the containers 
from above before experiments commenced. Images of 
the rods in motion were captured by a high speed dig- 
ital camera. The position and orientation of each rod 
was found by an image processing algorithm, developed 
in Matlab 6.5. The counting accuracy was close to 100% 
for low density systems and greater than 75% for samples 
as dense a.s (j) = 0.6. 

Unlike the infinite bulk modeled by Onsager's mean 
field theory rigid boundaries confine our experimen- 
tal system. We first investigate the FN transition away 
from these strong wall influences. By measuring the an- 
gular correlation of the rods with respect to the nearest 
wall [wall-rod correlation, 52 (j')], we distinguish rods af- 
fected by the side boundaries from those that were not. 
Briefly, we define an effective "bulk" when g2{r) ~ 0, 
which occurs at distances r/L > 0.5 from the wall. We 
will later discuss in detail 52 (^) and other parameters in 
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relation to ordering along the boundaries. 

To measure rod ordering in the bulk, we characterize 
alignment via the order parameter S = ( cos(20j) ), where 
6i is the angle of the «th rod with the nematic director. 
In practice, the order parameter and director are found 
for an arbitrary Cartesian coordinate system via the ma- 
trix Q — { 2ua{i)up{i) — Sap ), where Ua{i) and Ufj{i) are 
the coordinates of the unit vector specifying the direction 
of the ith rod. The eigenvalues of Q give ±S, and the 
eigenvectors give the director and its perpendicular vec- 
tor. A perfect nematic corresponds to S* = 1. Standard 
convention denotes an isotropic system as 5 = 0, but by 
the way S is defined, an isotropic system will have 5 > 
that approaches only as iV ^ (X) |9j . To distinguish rod 
ordering from this finite- number effect, we evaluated S 
of finite numbers of non-interacting rods using 2D Monte 
Carlo (MC) simulations, where the number of rods in the 
simulations matched those in experiments. 

We first tested the dependence of orientational order 
on area fraction, 0. Typical experimental results are il- 
lustrated in Fig. n where L/D = 40 in the R — 7 cm 
container. In the top of Fig. ^ the rods in the bulk 
(highlighted) appear to align as (p increases. The bottom 
of Fig. n shows S* as a function of as well as S from 
the MC simulations. For small (p, S is relatively low 
and coincides with S from simulations, indicating that 
the experimental system is isotropic. As 4> increases, S 
rises and deviates from the values derived from simu- 
lations, denoting the transition to an ordered nematic 
state. Conforming to the convention that S* = in an in- 
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FIG. 1: Dependence of bulk orientational order on rod density for a 
single rod aspect ratio L/D = 40 in the R = 7cm container. TOP: 
Representative images, for area fraction = 0.15 and 0.45. Dashed 
lines mark where 32(1") ~ 0, and rods with midpoints within these 
dashed lines were labeled "bulk" . BOTTOM: Order parameter, S, as a 
function of (p for experimental data (O)- Monte Carlo simulations (A) 
of non-interacting isotropically oriented rods are included to distinguish 
rod ordering from the finite-number effect on S. 



finite isotropic system, we henceforth subtract the results 
of the MC simulations from experimental data. 

We then fixed (p and tested the dependence of orien- 
tational order on rod aspect ratio, L/D. Experimental 
results for L/D = 20 and 60 at = 0.6 in the i? = 15 cm 
container are shown in Fig. EJtop). Though the same 
mass is present in both containers, the L/D = 60 rods 
appear more ordered than the L/D = 20 rods for this 
value of (p. The inset in the bottom of Fig. 13 shows S vs 
(p for all L/D and containers tested. As L/D increases, 
the transition to an ordered state shifts to lower (p. 

By plotting S versus the scaled area fraction, p — 
(pL/D, the transition curves for different L/D collapse 
to a single curve [Fig.jSJbottom)]. This scaling is consis- 
tent with Onsager's mean field result which assumes that 
the I-N transition stems only from steric interactions be- 
tween neighboring rods 0. The dashed line is a linear 
least square fitting of the data in the transition region. 
By extrapolating this line to 5 = 0, we estimate that 
the critical transition occurs at p = between 5 and 8. 
For comparison, Onsager's mean field model applied to 
an infinite 2D system of rods, predicts an I-N transition 
at p* = 37r/2 « 4.71 [l^, while 2D computer simulations 
with periodic boundary conditions estimate p* between 
6.5 and 8 We performed a similar analysis on our 
data that included both bulk and boundary layer rods. 
The measured S for the nematic state is systematically 
lower (due to cancellation by the alignment of the rods 
along orthogonal walls), but the location of p* remains 
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FIG. 2: Dependence of bulk orientational order on rod aspect ratio. 
TOP: Representative images, for L/D = 20 and 60 at </< = 0.3 in 
the R = 15 cm container. As in Fig. [3 rods within the dashed lines 
delineate bulk. BOTTOM: Consolidated plot of S vs. scaled area 
fraction, p = <j>L/D, for R = 7cm container: L/D = 20(A), 40{n); 
R = 15 cm container: L/D = 20(>), 60(O)l square container: 
L/D = 20(v), 40(0). Dashed line is a linear least squares fitting 
of the data in the transition region. By extrapolating this line to 
5 = 0, we estimate the transition density between 5 and 8. INSET: 
Same data plotted as S vs. 0. All solid lines are guides for the eye. 
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unchanged. 

While the side boundaries do not affect p*, the confin- 
ing walls profoundly influence rod ordering in other ways. 
Rods near a wall tend to align parallel along it. We as- 
sessed this alignment using the wall-rod angular correla- 
tion parameter g2{r) = (cos[2(-!5u) — ''^{r))] ), where d{r) 
is the angle of a rod located at distance r away from the 
wall, as described by the shortest line from rod midpoint 
to wall, and denotes the vector locally tangent to the 
wall. 

The top of Fig. El shows g2{r) for L/D = 20 in the 
R — 7 cm container at p = 1.5 (circles) and p — 9 (tri- 
angles). Experiments show an initially high g2(r) that 
decays as r increases. For comparison, the dashed line 
represents the MC calculated g2{r) for a single rod inter- 
acting only via hard core repulsions with the curved wall. 
Steric interactions with the wall significantly restrict the 
possible orientations available for a rod closer than L/2 
to the wall. Experimental values of g2{r) closely follow 
the single rod curve for small r; however, these values 
deviate as r increases. Furthermore, as p increases, g2i'f) 
diverges earlier and remains elevated for larger r. We also 
measured the number density, x(r), normalized with re- 
spect to bulk density, of rod midpoints (Fig. |2| bottom) . 
The measured x(r) is less than 1 for r/L < 0.5, which 
signifies a depletion of rod midpoints next to the wall. 
Below r/L ^ 0.4, x(^) increases and reaches higher val- 
ues with larger p. Finally, x(7') obtains a local maximum 
around r/L 0.6 before returning to bulk values. 

Onsager's approximation for excluded volume interac- 
tions has been applied to hard spherocylinders in contact 
with a hard flat wall by Poniewierski , and results have 
been confirmed by simulation Although this anal- 

ysis was performed for a 3D system contacting a planar 
wall, the experimental features in Fig. |21 show striking 
similarities. The profiles for g2{T) and x{f)^ however, 
can be difficult to interpret since geometrical constraints 
tightly couple rod orientation and position for rod mid- 
points close to the wall. By analyzing distances r' based 
on the closest end of each rod to a wall, rod orientation 
and position are decoupled, thereby simplifying g2{r') 
and profiles j^. Fig. Ol^top inset) shows a clear p 

dependent increase in (72 (j"') for rod ends near the wall. 
And, in contrast to the depletion layer seen for rod mid- 
points, a strong wetting layer becomes apparent, with an 
approximate four fold increase in x(^') for small r' which 
decays rapidly to bulk values around r'/L « 0.3 (Fig. 13 
bottom inset). 

Although many similarities exist between the hard core 
interactions of rods with straight versus curved bound- 
aries, wall geometry significantly affects 32 {r) . Rod align- 
ment with the wall does not penetrate as deeply into the 
sample center for curved boundaries and, interestingly, 
the thickness of this layer even regresses at very high p 
(data not shown). This seemingly unusual behavior can 
be understood by noting that for p ^ p* the system must 



accommodate both the orientation of the nematic direc- 
tor as well as the preferred orientation parallel to the 
surrounding boundaries. This frustration creates com- 
petition between the rods in the bulk with those at the 
boundary, which produces complex patterns, Fig.^fleft). 

Even though the local variations in the experimentally 
observed patterns are often on the order of rod dimen- 
sions, we find that these patterns are well approximated 
by minimizing a continuum free energy functional that 
describes elastic deformations in ordered liquid crystals 

0, ■ Consider a total free energy composed of bulk 
and wall contributions, T — Tbidk + J'waii- In the bulk, 
splay and bend deformations contribute to the distortion 
energy: 

Hulk = \ j Ki{\/ ■ nf + Ksin X V X n)^ dxdy, (1) 

where n{x,y) — [cosi?, sin?9] is the unit vector describ- 
ing the local nematic director, and Ki, are constants 
for splay and bend, respectively. Note that contributions 
from twist deformations are absent in this 2D system. 
The elastic constants relate to the molecular details of 
the nematic, generally becoming larger with increasing 
L/D or j^J. Often, however, the two constants are 
comparable in magnitude, and here we shall use the sim- 
plifying approximation, Ki = = JsT 0. To account 
for interactions favoring rod orientations locally parallel 
to the wall, we consider the line integral over the con- 
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FIG. 3: Wall-rod correlation function, g2('"), and normalized den- 
sity profile, x(^')' versus distance r from wall to rod midpoints for 
L/D = 20 in the R = 7cm container at p = 1.5 (circles) and p = 9 
(triangles). Data was smoothed by a moving average filter with a span 
of 3. Dashed line represents calculated g2('') for a single rod inter- 
acting with the wall via hard core repulsions. Insets show §2(1'') and 
X(r') versus distance r' from wall to rod ends. 
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tainer's perimeter P: 

T^aii = J^sm\^^ - i9) ds, (2) 

where C is a positive constant with units of Hne tension, 
and -dw describes the vector locally tangent to the wall, 
ritu . From Eq. |21 it becomes clear that J-waii depends not 
only on anchoring strength, C, but also on wall curvature 
and length. 

Fig.^Jright) shows the results of numeric minimization 
of the elastic free energy functional J- using a relaxation 
method with simulated annealing where the local 
tangent to each curve represents the orientation of n at 
that point. Similar patterns found in experiment (left) 
and calculations are grouped and arranged according to 
the relative ratio J-'buik/^waii- As nematic "stiffness" in- 
creases compared to wall anchoring {increasing K/C), 
^buikl^waii decreases, because the system reduces the 
cost of nematic deformations in the bulk at the expense 
of misalignment with the wall. 

Experimentally, by holding wall geometry constant 
and varying (j) and L/D, we changed nematic stiffness 
K and the ratio K/C, thus modifying J-buik/ ^waii 
Differing combinations of (f> and L/ D give rise to equiv- 
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FIG. 4: Elastic deformations of the nematic observed for p above 
the l-N transition. Experimental images in square and R = 15 cm 
containers (left) and results from numeric minimization of the elastic 
free energy functional, T, (right) are arranged with respect to relative 
bulk and wall free energies, Tbulk / ^wall ■ with Tbulk/^wall > 1, (a), 
~ 1, (b,c), and <C 1, (d). This ratio is modified experimentally by 
changing L/D or <)> (see text). A transition can be made from the 
patterns seen in (b,c) to those in (d) by either increasing L/D,(b), or 
<^,{c). 



alent patterns, as seen in Fig. ^b,c), corresponding to 
^buik/^waii ~ 1- By raising nematic stiffness though in- 
creasing L/D [Fig. 0fb)] or (j) [Fig- Etc)] i we transition 
to the pattern seen in Fig. Efd), where Tbuikj J'waii ^ 1- 
At the other extreme, where Tbuikj ^waii ^ 1, alignment 
with the wall is favored at the expense of nematic defor- 
mations. By changing wall curvature, we can strengthen 
wall anchoring and the penalty for wall misalignment, 
as seen in Fig. Efa). For a wide range of parameters, 
we not only find that the experimental patterns match 
the theoretical model results, but also that Tbuikl J'waii 
is sufficient to describe the whole range of nematic de- 
formations witnessed in experiments. However, as may 
be expected for confined nematics, we have also experi- 
mentally observed line and point defects in these phases. 

To conclude, we have shown that the macroscopic gran- 
ular rods in this steady state system share many proper- 
ties with thermally equilibrated lyotropic liquid crystals, 
such as an I-N transition and steric boundary effects. 
This system may also be useful for studying defects and 
other aspects of lyotropic liquid crystals that might be 
difficult to directly probe experimentally. 
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